First-principles study of phase stability of Gd-doped EuO and EuS 



O 

(N 

O 
CD 

Q 



CZ3 



J. M. AV'B S. V. Barabash, 2 V. Ozolins, 2 M. van Schilfgaarde, 3 and K. D. Belashchenko 1 

1 Department of Physics and Astronomy, Nebraska Center for Materials and Nanoscience, 
University of Nebraska, Lincoln, NE 68588, USA 
2 Department of Materials Science and Engineering, 
University of California, Los Angeles, CA 90095-1595, USA 
d Arizona State University, Tempe, Arizona 85284, USA 
(Dated: December 21, 2010) 

Phase diagrams of isoelectronic Eui-^Gd^O and Eui-^Gd^S quasi-binary alloy systems are con- 
structed using first-principles calculations combined with the standard cluster expansion approach 
and Monte-Carlo simulations. The oxide system has a wide miscibility gap on the Gd-rich side but 
forms ordered compounds on the Eu-rich side, exhibiting a deep asymmetric convex hull in the for- 
mation enthalpy diagram. The sulfide system has no stable compounds. The large difference in the 
formation enthalpies of the oxide and sulfide compounds is due to the contribution of local lattice 
relaxation, which is sensitive to the anion size. The solubility of Gd in both EuO and EuS is in 
the range of 10-20% at room temperature and quickly increases at higher temperatures, indicating 
that highly doped disordered solid solutions can be produced without the precipitation of secondary 
phases. We also predict that rocksalt GdO can be stabilized under appropriate experimental condi- 
tions. 
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I. INTRODUCTION 



The unique properties of gadolinium-doped europium 
chalcogenides make them attractive for spintronic and 
neutron detection applications. Doped EuO undergoes 
a spectacular metal-insulator transition near its Curie 
temperature T c , which is accompanied by huge magne- 
toresistance. In some samples, the resistivity changes 
by up to 13 orders of magnitude upon changing the 
temperature^ - — or by up to 6 orders of magnitude upon 
the application of an external magnetic field;^ Doped 
EuO can be epitaxially grown on Si and GaN sub- 
strates, and it demonstrates a very high spin polarization 
of the conducting electrons in the ferromagnetic (FM) 
state. 4 This half-metallic behavior suggests applications 
of doped EuO as a spin injector material in spintronic 
heterostructures. On the other hand, Gd-doped semi- 
conductors are appealing as neutron-absorbing materi- 
als for solid-state neutron detection technology due to a 
very high neutron absorption cross-section of the 157 Gd 
isotope^— 

The potential applications of Gd-doped Eu chalco- 
genides depend on their phase stability. No phase dia- 
grams are available for Eu-Gd oxides; the data for Eu-Gd 
sulfides is limited to the EuS-Gd2S3 isoplethal section. 9 
The experimental studies of pure and Gd-doped EuO 
and EuS give only indirect information about the phase 
stability while focusing on other properties, which we 
briefly review here. Pure EuO is a rare FM insulator 
with a rocksalt structure, whose optical absorption gap 
increases monotonically from 0.95 eV at K to 1.12 eV at 
300 KiS. A divalent Eu ion in EuO has the s S 7/2 ground 
state configuration and the magnetic moment of 7/j.b due 
to the half-filled 4/ shell. The Curie temperature rises 
sharply with Gd doping from T c = 69 K in pure EuO up 



to T c = 170 K at optimal doping of about 4% 
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miscibility of Gd in EuO is unknown but is expected to be 
finite, particularly because GdO has not been observed 
in the rocksalt structure. In fact, the common Gd oxide 
has Gd203 stoichiometry (crystallizing in three different 
phases^), but a tentative observation of zincblende GdO 
has also been reported.— On the other hand, both EuS 
and GdS are stable in the rocksalt structure, and they 
can form a continuous range of rocksalt solid solutions at 
all concentrations^ EuS is an insulator with the absorp- 
tion gap of 1.64 eV, which is FM below T c = 19 GdS 
is an AFM metal with a Neel temperature of 58 K.— 

Several recent first-principles studies^— have focused 
on the electronic structure, magnetic interaction, and 
other properties of pure Eu monochalcogenides. The 
interplay between the magnetic ordering, spectral, and 
transport properties of doped EuO was also studied using 
model Hamiltonians.— ~— The structural phase stability 
of the (Eu,Gd)0 and (Eu,Gd)S quasi-binary alloys thus 
remains unexplored. 

The purpose of the present study is to give a theoreti- 
cal perspective of the phase identity, stability, miscibility 
and other properties of Gd-doped EuO and EuS along 
the EuO-GdO and EuS-GdS isoplethal sections. Based 
on the comparison with the self-consistent quasiparticle 
GW calculations and with available experimental data, 
we adopt the generalized gradient approximation (GGA) 
with the Hubbard U corrections applied only to the rare- 
earth 4/ orbitals. We then apply the first-principles 
cluster expansion technique and subsequent Monte Carlo 
simulations to construct the phase diagrams of (Eu,Gd)0 
and (Eu,Gd)S quasi-binary alloys. We find that, de- 
spite the isovalency of the two alloy systems, the result- 
ing phase diagrams are quite different. In particular, we 
predict that the oxide system has two yet-unobserved or- 
dered phases with 1:1 and 1:2 Gd-to-Eu ratios, which be- 
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come thermodynamically stable below ~ 900 and ~ 500 
K; moreover, we predict that rocksalt GdO can be sta- 
bilized in a narrow range of oxygen pressures. On the 
other hand, such 1:1 and 1:2 phases do not appear in the 
(Eu,Gd)S system, and moreover a different ordered phase 
with a 2:1 Gd:Eu ratio is very near the tie-line of the end 
compounds EuS and GdS. We further analyze the role of 
the chemical and deformation-mediated interactions and 
find that the qualitative difference between the oxide and 
sulfide systems is mainly due to the contribution of local 
anion relaxations. 

The paper is organized as follows. The methodolog- 
ical and computational details are described in Section 
HTl Section Mil presents the calculations of band struc- 
ture and elastic properties of end compounds. Compari- 
son with experiment and with GW calculations serves to 
justify the adopted GGA+C7 approach. Section [TV] dis- 
cusses the magnetic ordering energies and their relevance 
to phase stability. The configurational Hamiltonians are 
described in Section |V] followed by the calculations of 
phase diagrams in Section IVI1 The anion-mediated de- 
formational interaction mechanism is discussed in Section 
IVII1 and finally we summarize in Section IVHIl Some 
technical details, including the extraction of paramag- 
netic formation enthalpies and the structural information 
for the predicted compounds, are included in the Appen- 
dices. 



II. COMPUTATIONAL APPROACH 
A. Total energy calculations 

Total energy calculations for all ordered compounds in 
this study were performed using the projected augmented 
wave (PAW) metho d 26 ' 27 and generalized gradient ap- 
proximation (GGA) of Perdew-Burke-Ernzerhof^ with 
the Hubbard U correction 2 ^ for the 4/ orbitals, as imple- 
mented in the VASP packaged An energy cutoff of 500 
eV was used for the plane-wave expansion of wave func- 
tions, and total energies were converged to within a few 
meV per atom with the density of the k-point mesh no 
lower than 0.01 A -3 including the V point for Brillouin- 
zone sampling. 

From the energetics of the individual atomic levels, it 
is clear that in both oxides and sulfides, the 4/ bands are 
half-filled, whether originating from Eu or Gd, whereas 
the binding energy of the Gd 4/ states is much larger 
compared to divalent Eu. Electron doping through the 
addition of Gd fills the conduction-band states, which 
are spin-split by the exchange interaction with the 4/ 
shell. The shallow core 5s 2 and 5p 6 states on both Eu 
and Gd are included in the valence basis set. The half- 
filled and strongly-localized 4/ orbitals in both Eu and 
Gd were treated within the GGA+U approach.— To cal- 
culate the value of J = 0.6 eV, we used the constrained 
occupation method 3 -' 32 by considering the 4/ states as 
an open-core shell and finding the total energy difference 



between the 4/j 5 and 4/.j? 4/j states. This calcula- 
tion was performed using the full-potential linear aug- 
mented plane wave (FLAPW) method implemented in 
the FLEUR packaged The value U = 5.3 eV calculated 
in a similar way, when used in the GGA+U calculation, 
leads to the 4/ states being too shallow with respect to 
the conduction band of EuO. This discrepancy is due to 
the underestimation by GGA of the intrinsic insulating 
gap between the 0-2p and Eu-5d states. This underes- 
timation (which is not sensitive to U) brings the 0-2p 
states too close to the Eu 4/ states; their hybridization- 
induced repulsion pushes the 4/ states up and reduces the 
band gap. Therefore, we adjusted the value of U empir- 
ically and found that U = 7.5 eV results in good agree- 
ment with optical absorptio n 34 ' 35 and photoemission^ 6 . 
measurements, as illustrated in section Mil The Hub- 
bard U corrections are not used for the 5d orbitals on Eu 
and Gd, as justified below in sections Til Bl and Hill 

The proper treatment of the chemical disorder in many 
/ and d electron systems is notoriously difficult, because 
of the dependence of the structural energy on the / (or 
d) orbital orientation, and due to the large spurious self- 
inter acion present in the popular GGA+U versions for / 
(or non-time-reversible d) orbitals^ Fortunately, this is 
not an issue in (Eu,Gd)0 and (Eu,Gd)S, since the half- 
filled / shells do not exhibit such an orientational energy 
dependence. Indeed, for all considered cases, we found 
the lowest energy electronic configurations to correspond 
to exactly seven co-aligned / electrons on each cation 
atom, forming a rotationally-invariant shell. At the same 
time, it was not unusual for our early test calculations 
to lead to electronic states with different /-electron oc- 
cupation, with energies higher by 0.15. . .3.0 eV/cation. 
Some of such higher-energy states were clearly identifi- 
able as having a Gd electronic configuration of 4f 8 5d° 
instead of 4/ 7 5d 1 . We found that the appearance of such 
states was either a failure of the residual minimization 
method (RMM),— or an artifact of an insufficient initial 
/-electron spin polarization. By using the Davidson min- 
imization algorithm 3 - and by assigning the initial on-site 
spin polarization of 10 \xb per cation (which is partially 
assigned to the d and s electrons, resulting in the desired 
occupation of 7 co-aligned / electrons), we could avoid 
such high-energy electronic configurations. The only re- 
maining degree of freedom for the final electronic state 
corresponds to the magnetic ordering of the fully polar- 
ized cations, which is well controlled by the signs of the 
initial on-site spin polarizations and is further discussed 
in Section HVl 

The spin-orbit corrections are often appreciable for 
the band energies of heavy-element compounds at high- 
symmetry k-points, but are typically minor for the forma- 
tion energies, which include contributions averaged over 
different bands in the entire Brillouin zone. Indeed, we 
checked the effects of the spin-orbit coupling on the for- 
mation enthalpies AiJoGA for the key ordered structures 
that we identify below in Sec lVIl as the stable ground 
state compounds. We found that the AHqqa change 
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by at most 2 meV/cation (see Table IIIII below). We 
therefore performed all our routine calculations in the 
scalar-relativistic approximation, without including the 
spin-orbit coupling. 



B. Benchmark GW calculations 

It has been suggested by some author o 19 ' 20 that the 
addition of Hubbard corrections for the empty 5d states 
may be necessary for the correct description of the con- 
duction band, and in particular for the determination 
of the character of the band gap (direct or indirect) of 
europium chalcogcnides. However, these corrections are 
arbitrary unless a reliable benchmark is used to select 
the Hubbard parameters. We resolve this issue by calcu- 
lating the band structure of EuO using the quasiparticle 
self-consistent GW (QS GW) approximation. Here, we 
discuss the QSGW methodology, and later in Section Hill 
we use QSGW to show that the conduction band struc- 
ture comes out almost exactly right in the GGA+17 cal- 
culation with U applied only to the 4/ orbitals. 

QSGW has been shown to be a reliable predictor of 
materials properties for a wide range of compounds com- 
posed of elements throughout the periodic table42r— 
Nevertheless, prior experience has revealed certain kinds 
of systematic errors. The correction of these errors makes 
minor adjustments to weakly correlated materials sys- 
tems, and somewhat stronger adjustments for more cor- 
related materials. There are two highly systematic errors 
that affect the band structure of EuO. 

First, bandgaps in semiconductors such as GaAs, and 
insulators such as SrTiC>3 are systematically overesti- 
mated a little. The same effect is seen in the spd subsys- 
tem of EuO. Second, shallow core-like levels, such as the 
highest occupied d levels in Zn,Cd,Cu,Ag,Au, and so on, 
are systematically predicted to be too close to the Fermi 
level, typically by < 0.5 eV. This error is seen in the / 
subsystem of EuO, as we will discuss. 

Both types of errors are highly systematic in sp and d 
systems, and discrepancies with experiment in 4/ com- 
pounds are consistent with these errors4^ To a large ex- 
tent, the first error can be simply explained through the 
random phase approximation (RPA) to the screened in- 
teraction W, which can be understood as follows. The 
RPA bubble diagrams do not include electron-hole inter- 
actions in intermediate states in the calculation of the 
irreducible polarizability n(q, w) and thus the dielectric 
function e(q, oS). Short-range attractive (electron-hole) 
interactions induce the red shifts in Ime(q,w) at ener- 
gies well above the fundamental bandgap; see e. g. Fig. 6 
in Ref.HcJ Ladder diagrams are sufficient to remedy most 
of the important errors in II(q,u>), as was demonstrated 
rather convincingly in Cu204^ 

Inclusion of these contributions increases the static 
dielectric constant £„, as can be readily seen through 
the Kramers-Kronig formula relating the real and imag- 
inary parts of e. Remarkably, e^ calculated by the 



RPA in QSGW is underestimated by a nearly uni- 
versal factor of 0.8, for many kinds of insulators and 
semiconductors^ 2 - including transition metal oxides such 
as C^O^S, SrTiOg^, Ce02, and sp semiconductors^ 
Because e is systematically underestimated, W = e~ l v 
(where v is the Coulomb interaction), the self-energy 
S = —iGW, and the quasiparticle excitation energies 
are systematically overestimated. 

The second kind of error cannot be explained in this 
way. QSGW pushes down the semicore d level in Zn 
(or / level in Eu) rather strongly relative to the LDA; 
however, the shift is slightly too small to agree with ex- 
periment. As we have noted, W calculated by QSGW 
is already too large: reducing W reduces this correction. 
This implies that the error should be attributed to the 
other approximation in GW theory, namely the omission 
of the vertex T in the formally exact self-energy, GWT. 

Both kinds of corrections to the GW approximation 
are difficult to carry out in practice. As for the correction 
to n, we have noted that simply scaling E by 0.8 largely 
undoes this error, in a wide range of systems. We make 
such a scaling here, to correct the spd subsystem. 

Whether or not E is scaled, the Eu d-f gap is too small: 
when E is scaled it comes out approximately zero, in 
contrast to the observed gap of about 1 eV. It is expected 
that the vertex in GWT will largely just induce a shift 
in the Eu / state. Anticipating this, we included an ad 
hoc addition to the QSGW potential for EuO, essentially 
doing a QSGW+U calculation with U = 0.816 eV. The 
value of U is adjusted to make the d-f gap coincide with 
the experimental number. With these corrections, we 
anticipate the QSGW method to yield highly accurate 
band structures, discussed below in Sec. IIIII 



C. Cluster expansion and Monte Carlo simulations 

To identify the thermodynamically stable phases and 
their range of stability, the standard cluster expansion 
(CE) formalis m 48 i 49 coupled with the ground state search 
and Monte Carlo simulations were employed, using the 
routines implemented in the Alloy- Theoretic Automated 
Toolkit (ATAT) packagei 50 i 51 We consider quasi-binary 
substitutional alloys, assuming that the anion sublattice 
is fully occupied by the chalcogenide atoms of one kind 
(either O or S). Throughout the paper, we use the terms 
"structure" and "ordering" to refer specifically to the or- 
dering of Eu and Gd cations within their own (fee) sub- 
lattice, unless indicated otherwise. Our CEs cover FM 
compounds, with other cases considered separately, as 
detailed in Sec. IIVI The cell size and shape, as well as all 
the atomic positions, were fully relaxed for all structures 
using the conjugate gradient algorithm, starting from the 
ideal rocksalt lattice. In view of prohibitive computa- 
tional cost, we did not consider the phonon contribution 
to the formation enthalpies, which can somewhat modify 
the phase diagrams at elevated temperatures. 

The relaxed formation enthalpies are parameterized by 
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a configurational cluster-expansion Hamiltonian AHce'- 

AH C E(<r) = J2 J f D f n W ' C 1 ) 

/ 

where the occupational degrees of freedom are described 
by a configurational vector er (a particular decoration 
of the cation sites of the rocksalt lattice by Gd and Eu 
atoms), Jf is called the effective cluster interaction (ECI) 
for a cluster figure / with Dj as the figure's symme- 
try degeneracy per site, and n(<x) is the configuration- 
dependent correlation function in the interaction clus- 
ter. In practice, a finite number of terms nf is kept in 
the expansion ([T]) , and the expansion becomes exact*^ as 
Uf — > oo. 

The ECI values J/ are determined by fitting to a set of 
N- m "input" formation enthalpies AHqca (""in)- Mn was 
iteratively increased by performing GGA+U calculations 
for new structures <7 ln based on the CE predictions, until 
a desired CE accuracy was reached, in particular estab- 
lishing an agreement between the ground state predic- 
tions of the final CE and GGA+J7. (For the oxide system 
such full consistency was established only within a target 
concentration range discussed in Sec I VI Al ) In order to 
evaluate the predictive power of the cluster expansion, 
the "leave-one-out" cross-validation (CV) score was cal- 
culated using the procedure implemented in ATAT^S for 
each (To out of the N- m input structures, a separate fit- 
ting of the Jf values was performed with that one struc- 
ture excluded from the fitted set. The actual energy 
Ai?GGA(co) of the excluded structure was then com- 
pared with the prediction AH^ a " (<x ) of this "leave-one- 
out" fitting, and the difference was averaged over all the 
N[ n choices of <tq . Unlike the conventional mean-squares 
fit error, which monotonously decreases upon increasing 
n/, the CV score is designed to measure the predictive 
power of the cluster expansion and has a minimum for a 
finite rif value, diverging if rif becomes too large. 

Once the ground states for a range of concentrations 
have been identified, the ones that are stable at T = 
K are determined by the convexity condition; the given 
structure at concentration x is stable if it lies below any 
straight line connecting other compounds at concentra- 
tions X\ < x and X2 > x. (The convex hull of the set of 
points in the AH vs x plot represents the full range of 
enthalpies that may be achieved by the system. Stable 
compounds are those that form the vertices on the lower 
boundary of this convex hull.) For each stable compound 
one can define its "energetic depth" 6, i. e. the amount 
by which its energy would increase if it were decomposed 
into two stable compounds that are closest to it in con- 
centration. 

The phase diagrams were computed using semi-grand 
canonical Monte Carlo (MC) simulations (i. e. with the 
varying number of Eu and Gd cations) and the Metropo- 
lis algorithm implemented in ATAT. For the final phase 
diagram construction, we used an 18 x 18 x 18 super- 
cell based on primitive fee translations (5832 cations in 
the simulation box). In the case of (Eu,Gd)0, we esti- 



mated the effects of the finite size and the commensu- 
rability with other ground state structures by also using 
12 x 12 x 12 and 15 x 15 x 15 simulation boxes. We found 
that the 18 x 18 x 18 box was commensurate with all 
the ground-state structures that have ordering temper- 
ature T or d above 400 K, although this did not hold for 
some structures stable at lower T (either identified by 
the CE ground state search or directly observed in MC 
simulations). All such low-T or( j structures in the oxides 
are limited to the Eu-rich composition range indicated 
below in the phase diagram. A series of simulations was 
performed at fixed chemical potentials with temperatures 
varying in 2-5 K increments. The equilibration and sam- 
pling passes were done with 1000-5000 Monte Carlo steps 
(flip attempts per site); longer runs of 10000 steps were 
performed for particularly difficult regions. The phase 
boundaries were then found by identifying the disconti- 
nuities (or cusps) in the dependence of the average con- 
centration and enthalpy on temperature. 



III. BINARY CHALCOGENIDES: PROPERTIES 
AND ELECTRONIC STRUCTURE 

Figure [1] shows the electronic band structures of EuO 
and GdO calculated using both GGA+U and QSGW ap- 
proaches. The densities of states of all four end com- 
pounds calculated in GGA+C/ are shown in Fig. [2j 

When Hubbard U is applied only to the rare-earth 4/ 
states, the conduction band structure of both EuO and 
GdO is in very good agreement between GGA+U and 
QSG1U, the band gap being indirect. On the other hand, 
the addition of the Hubbard U correction to the Eu 5d 
orbitals proposed in Ref. [H pushes the bands up at the 
X point but not at T, resulting in a direct band gap 
in disagreement with QSGW. Therefore, as mentioned 
above, we use GGA+C/ with U applied only to the rare- 
earth 4/ states. Physically, this is reasonable because 
the 5d states of either Gd or Eu are rather delocalized 
and have a sizable bandwidth, being thus amenable to 
treatment within GGA. 

The main difference between the GGA+C/ and QSG1U 
band structures is the position of the oxygen p states 
in EuO. In GGA+Z7, they lie approximately 2.5 eV 
too high, resulting in a stronger hybridization of the 
majority-spin O states with the occupied 4/ states. This 
hybridization leads to a repulsion of these states from 
each other, and to a large spin splitting of the oxygen 
p states. The QSGW results are in very good agree- 
ment with photoemission measurements^ By contrast, 
for GdO the oxygen p states come out only a little too 
shallow. This is reasonable, because the gap between 
the oxygen p and rare earth 5d states in GdO is not ac- 
companied by a discontinuity of occupation numbers; the 
latter is responsible for the band gap problem in semi- 
conductors. We also note that QSGW significantly over- 
estimates the energy of the unoccupied 4/ states, which 
appears to be its universal feature^. This error is imma- 
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TABLE I: Calculated (this work, marked "theory") and ex- 
perimental (when available, marked "exper.") lattice parame- 
ters, elastic constants, and the bulk modulus B of the binary 
rocksalt compounds. All elastic constants are given in Mbar 
units. 



Compound 


a, A 


Cu 


Cl2 


C44 


B 


EuO the ° ry 
exper. 


5.18 
5.14 


1.89 
1.9" 


0.62 

0.42(8)" 


0.78 
0.54" 


1.04 
0.92(6);" 1.10 6 


GdO theory 


4.92 


3.54 


0.63 


0.72 


1.60 


EuS the ° ry 
exper. 


6.02 
5.97 


1.47 
1.3" 


0.24 
0.11(8)" 


0.39 
0.27" 


0.65 

0.51;" 0.61; 6 0.72 c 


GdS theory 
exper. 


5.62 
5.56 


3.06 


0.33 


0.36 


1.24 
1.20 



a Ref. HI 
b Kei. |H 
c Ref. |H 



terial for our purposes. 

The choice of U = 7.5 eV and J = 0.6 eV simulta- 
neously produces the splitting between the occupied and 
unoccupied 4/ states of U + 6 J = 11.1 eV in good agree- 
ment with photoemission and inverse photoemission mea- 
surements for the Eu metal^ the optical band gap at the 
X point of 0.94 eV in EuO consistent with the value of 
0.95 eV measured at zero temperature; 34 ' 35 and the equi- 
librium lattice constant a = 5.182 A in good agreement 
with the experimental value of a = 5.144 A. The appli- 
cation of the same U and J values to EuS leads to the 
optical gap at the X point of 1.52 eV consistent with the 
zero-temperature value of 1.51 eVJ£ Due to the similar 
nature of the half-filled 4/ orbitals, these semiempirical 
U and J values were applied to both Eu and Gd 4/ states 
in the oxide and sulfide systems. Both GdO and GdS arc 
metallic, as expected. 

Table U includes the lattice parameters and elastic con- 
stants calculated in the FM state, along with the avail- 
able experimental data for EuO, EuS, and GdS. The 
lattice constants are slightly underestimated by approx- 
imately 1%. The bulk moduli agree within the un- 
certainty of the experimental data. The elastic con- 
stants obtained by ultrasonic measurements have only 
been reported by one groupi 5 - 3 - The calculated Cn con- 
stant agrees well with this measurement for both EuO 
and EuS. The C44 constant is overestimated by approxi- 
mately 45%. The C12 constant is also overestimated, but 
comparison is hindered by a very large experimental un- 
certainty. Note, however, that the bulk moduli obtained 
from these measured elastic constants are among the low- 
est ones reported in the literature; it is possible that these 
measurements are affected by off-stoichiometry. There- 
fore, it is unclear whether the disagreement in C44 and 
C12 is due to the inaccuracy of the GGA+U method or 
to experimental artifacts. 



IV. MAGNETIC ORDERING 

The magnetic ordering temperatures of pure EuO, EuS 
and GdS (69, 19, and 58 K, respectively), as well as for 
the entire range of solid solutions, are well below room 
temperature. Therefore, all structural phase equilibria 
involve paramagnetic (PM) phases with randomly ori- 
ented local moments on the Eu and Gd atoms. In princi- 
ple, the PM enthalpies for the input structures can be es- 
timated by fitting a number of magnetic configurations to 
a Heisenberg Hamiltonian and taking the constant term 
as the PM energy. However, doing this for more than a 
few simple ordered compounds is computationally pro- 
hibitive. Fortunately, relatively low magnetic ordering 
temperatures suggest that the use of ground state forma- 
tion enthalpies to study configurational thermodynamics 
should not lead to large errors. Nevertheless, in this sec- 
tion we perform a few checks and discuss the possible 
modifications introduced in the cluster expansion by the 
replacement of FM formation enthalpies by the PM ones. 

We considered several collinear magnetic configura- 
tions in seven simple (Eu,Gd)0 and four (Eu,Gd)S 
compounds;~ these results are summarized in Table [TTJ 
Like in the rest of our study, the scalar-relativistic ap- 
proximation was used here. Where a comparison can be 
made, we found good agreement with other published 
data- 1 ^ For the (Eu,Gd)0 system the FM state al- 
ways has the lowest energy, except in pure GdO where 
the AFM phase with the ordering vector along [111] is 
slightly lower (by 0.6 meV/cation) than the FM phase. 
This indicates that restricting our CE study to the FM 
compounds is sufficient to yield an accurate description 
of the thermodynamic phase stability at T = K. 

One can argue that FM enthalpies can also be used 
to predict the phase stability of oxides at higher tem- 
peratures. As mentioned in the introduction, the T c of 
Gd-doped EuO rises sharply to about 170 Katxw 0.04 
and then slowly decreases (see e.g. Ref. |l~2"[). However, 
the magnetization curves for these higher-T c alloys have 
a distinctive double-dome shape. Specifically, the mag- 
netization drops to a fairly small value at temperatures 
close to T c of pure Eu ("main dome"), and extends a rel- 
atively weak tail up to the actual elevated T c . While the 
mechanisms of this behavior are not completely under- 
stood, it is fair to assume that the dominant part of the 
PM-FM enthalpy difference is released in the tempera- 
ture range of the "main dome." Since the characteristic 
temperature of this feature does not strongly depend on 
the doping level, we expect that the PM-FM enthalpy dif- 
ference is a featureless function of the concentration, and 
that it likely does not exceed 10 meV^ This conclusion 
is generally consistent with the data in Table |Hj Such 
correction is not likely to lead to significant changes in 
the phase diagram, and we therefore use FM enthalpies 
for oxides in the following. 

In the sulfide system the situation is different for two 
reasons. First, as we discuss in Section IVI B[ the FM CE 
predicts an ordered EuGd2S3 structure with a high or- 
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(c) GGA+U : FM GdO (d) GW : FM GdO 




lt xw r lt xw r 



FIG. 1: (Color online) Electronic band structure of EuO calculated (a) in GGA+U and (b) in QSGW, and of ferromagnetic 
phase GdO, calculated by (c) GGA+U and (d) QSGW. Black (red) curves correspond to the majority (minority) spin bands. 
Energy is referenced either from the valence band maximum (VBM) or from the Fermi level (Ef). In GdO, states near Ef at 
F are of Gd d character; the band whose value is —1.5 eV at r is of s character. If correlations were strong, the d band would 
become narrow and possibly shift relative to the s band. As can be seen, GGA and QSGW are very similar for these bands: 
GGA and QSGW differ mainly in the positions of the O 2p bands, at around —6 eV. That O 2p states shift downward relative 
to GGA (or LDA) seems to be a universal property of oxide insulators. 



dering temperature, which is only marginally stable with 
respect to the pure FM EuS and GdS. Second, the mag- 
netic order changes from FM in EuO to AFM in GdS, 
and the latter has a relatively high FM-AFM energy dif- 
ference of 13.5 meV/cation (see Table [XT]) . Therefore, 
for EuS, EuGd 2 S 3 , and GdS we have estimated the PM- 
FM enthalpy differences, which are —5.3, +1.3, and —0.6 
meV/cation, respectively (see Appendix [SJ. These dif- 
ferences are sufficient to make the EuGd2S3 structure 



marginally unstable. The effect on the phase diagram is 
considered in Section IVIBI 



V. CLUSTER EXPANSIONS 

In this section, we characterize the CEs obtained at the 
end of the iterative CE construction procedure. These 
CEs are used in Sec lVII to evaluate the phase stability of 



> 

09 



00 

O 
Q 





1 


V ' 'euo'fm' 1 ' 1 : - 

\l , 1 r — — 1 — V—LA j 


1,1,1,1,1,1,1,1 




1 ^n*^-A 


i 


1 1 1 1 1 1 1 1 1 1 1 1 

\ EuS FM 




— ■ — I 1 1 L— ^i^_J^— ^- ~'L-~I 

- 1 , 1 I 1 , 1 , 1 I 1 , 1 , 1 1 





10 



10 



10 



12 




-2 

E-E P (eV) 



FIG. 2: Density of states (DOS) of FM EuO and EuS, AFM 
type II GdO, and AFM type II GdS calculated in GGA+U. 
For the AFM II phases, the solid line shows the total DOS in- 
cluding both the cation and the anion contributions, whereas 
the dashed lines show the partial majority and minority-spin 
DOS from the Gd cations. Energy is referenced from the 
Fermi level Ef- 



Eui-^GdrO and Eui_ x Gd x S rocksalt alloys. 



A. Eui-.Gd.O 

The initial input set for the self-consistent CE con- 
struction included N? n = 26 structures including all the 
atoms up to 4 cations per cell, except two such structures 
at EuGd30*4 composition. The final CE for Eui-^Gd^O 
has an input set of N- ln — 148 structures (identified 
throughout the CE iterations as potential ground states 
or otherwise as structures important for the CE accu- 
racy), and uses an ECI set of 8 pairs, 12 triplets and 
16 quadruplets. The predictive power of this CE is esti- 
mated by the CV score as 5.8 meV/cation, whereas the 
root-mean-square fit error for the input structures is only 
2.8 meV/cation. The ground state search was performed 
among all the structures up to 30 atoms per cell (~ 2 15 
configurations) . 

The ECIs values J/ for Eui_ x Gd x O as a function of 
the effective radius (the longest intersite distance in the 
cluster /) are shown in Figure |3] as the crosses connected 
by the red line. The leftmost panel in Fig. [3] displays the 
pairwise ECIs up to the eighth nearest-neighbor in the 
cation sublattice, and the right two panels correspond 
to triplet and quadruplet ECIs, respectively. By far the 
largest ECI is due to the second nearest-neighbor pair. 



TABLE II: Formation energies (in meV/cation) of select or- 
dered (Eu,Gd)0 rocksalt compounds with different magnetic 
orderings including FM, AFM with two different orientations 
of the layers of co-aligned spins, and ferrimagnetic [where ap- 
plies, with different arrangement of (111) layers of up and 
down spins relative to the majority (A) and minority (B) 
cation species]. Magnetic orderings requiring large cells were 
not considered, as indicated by dashes. 



Compound FM 



AFM 
(111) (001) 



ferrimagnetic (111) 
Af AfBJ, AtA|Bt 



EuO 
GdO 
EuS 
GdS 

Ll EuGd0 2 
Lli EuGd0 2 
Lli EuGdS 2 
C6 Eu 2 Gd0 3 
C6 EuGd 2 3 
C6 EuGd 2 S 3 










60.6 
-59.3 
8.9 
-41.5 
-32.9 
-0.4 



12.2 
-0.6 
0.1 
-13.5 

-47.9 
7.5 



15.8 
15.1 
3.4 
3.8 
78.1 



8.2 
1.6 



-36.6 
-21.4 
-3.6 



-33.5 
-17.0 
-4.6 



FM 



ferrimagnetic (001) 
AtA^AtB^ AtA|A|Bt AtA4.A4.B4, 



Zl Eu 3 Gd0 4 5.8 



17.1 



29.9 



22.4 



It is positive, and its magnitude (17.6 meV) is more than 
three times greater than that of the negative first nearest- 
neighbor ECI (—5.4 meV). Overall, the pair interactions 
are stronger than all other cluster ECIs. As explained 
in Section [VIII below, the large positive ECI for second- 
nearest neighbors is due to the significant energy gain 
from displacing O atoms towards Gd in the Eu-O-Gd 
double-bond patterns along any of the [001] directions. 
This ECI is the main driving force for the stabilization of 
the Lli cation ordering that we find below in Sec lVI Al as 
indeed could be expected^, from a simple CE involving 
only pair interactions J"^ r and Jp^™ between nearest- 
nighbor and second-nearest neighbor atoms of the same 
magnitude as in our actual CE. 



B. Em zGd^S 

The cluster expansion has been constructed with a CV 
score of 4.0 meV/cation using N- m = 49 input structures 
and an ECI set of 8 pairs, 7 triplets, and 3 quadruplets. 
The set of structures used for ground state search con- 
tained ~ 2 10 configurations, and the root-mean-square 
fit error for the input structures is only 2.2 meV/cation. 

The pair and many-body ECIs for the oxide and sul- 
fide systems are qualitatively similar (see Fig. [3]) . How- 
ever, we shall see in Sec lVIl that there is a striking dif- 
ference in the ground-state energetics, convex hull struc- 
tures, and the phase diagrams of the oxides and the sul- 
fides. This difference stems from the quantitative change 
in the ECIs: while the positive second nearest-neighbor 
pair ECI is still the strongest one in the sulfide system, 
it is reduced by about 40% compared to the oxides. This 
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FIG. 3: (Color online) Three different types of ECI parame- 
ters as a function of effective radius in units of the nearest- 
cation-neighbor distance r nn in the cluster expansions for 
Eu^Gdi-^O, denoted by red crosses, and EuzGdi-^S, de- 
noted by blue circles. 



reduction results in the destabilization of the Lli cation 
ordering, as we discuss in Section IVlII 



VI. PHASE STABILITY 

A. Eui-.Gd.O 

Figure [4ji shows the calculated formation enthalpy di- 
agram for the entire range of x in Eui-^Gd^O. Each 
green dot represents the composition and the CE forma- 
tion enthalpy AHce{g) of one of the 2 15 structures a 
used in the ground state search. The fitted AHcE^in) 
values for the 148 input structures are shown by the red 
crosses, and the input AHGGA{&in) are shown by open 
black circles. The input set contains eight stable com- 
pounds, serving as the breaking points on the GGA+U 
convex hull. (This convex hull is shown by the black 
solid line in Fig. 2J) These compounds represent tenta- 
tive ground state predictions, in a sense that they all were 
predicted as ground states throughout the CE iterations 
and confirmed by direct GGA+U calculations. Most of 
them have small energetic depths (see Table the 
phase diagram calculations show that only two of these, 
with x = 0.5 and x = 1/3, can appear at reasonably high 
temperatures. Therefore, it is not necessary to insist on 
the precise prediction of the ground state sequence at 
x < 1/3. Indeed, we found that the convex hull cor- 
responding to our final CE fitted Hamiltonian was not 
the same as that for the input set (i.e., it changes even 
throughout the final CE iterations). However, full self- 
consistency was achieved for x > 0.3, allowing us to iden- 
tify the compounds with x = 1/3 and 0.5 as unambigu- 
ously established ground states. The structural informa- 
tion for these two compounds is given in Appendix [B] 



Despite the slight disagreement regarding the identity of 
the low-temperature ground states at x < 0.3, the overall 
convex hull shape and the AHce values agreed well be- 
tween the input (GGA+U) and the final predicted (CE) 
convex hulls, and the identity, formation enthalpies and 
energetic depths of the x = 1/3 and 0.5 ground states 
are accurately reproduced. Note also that the resulting 
CE fitting does not have an exact meaning of a forma- 
tion enthalpy, because the vanishing of the fitted quan- 
tity for the end members is not enforced. In particular, 
pure EuO has a spurious "formation enthalpy" of —10.6 
meV/cation in the CE fitting. However, the shape of the 
ground-state convex hull suggests that the fit error for 
pure EuO should not affect the phase diagram signifi- 
cantly at x > 0.1 where the fitting is quite accurate. 



(a) Eu, Gd O 

* ' 1 -Y Y 



(b) Eu,_ Gd x S 
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FIG. 4: (Color online) Formation enthalpy per cation vs com- 
position x for all distinct cation orderings within rocksalt 
structure for (a) EuzGdi-zO with up to 30 atoms per unit cell 
and (b) EuzGdi-zS with up to 20 atoms per unit cell. The 
black open circles are the first principles inputs, red crosses 
are the fitted CE values for the input structures, and the green 
dots are the predicted AHce for all other structures. 

The calculated formation enthalpies for both the 
unambiguously established and tentatively predicted 
ground-state Eui-^Gd^O compounds identified through- 
out the CE iterations and confirmed by direct GGA+U 
calculations are listed in Table IIIII along with their CE- 
fitted values, as well as the energetic depths <5 found from 
both calculated and fitted formation enthalpies. The fit 
error for all these compounds is less than 4 meV/cation. 
The ground state with the largest 5 has Lli structure, 
which is an AiBi (111) superlattice, i.e. it is formed by 
alternating (111) layers of pure Eu and Gd. Among the 
six compounds with x < 1/3 there are three other pure 
Eu/Gd superlattices: Ai 4 B (134), Ai 2 B (124), and A 8 B 2 
(123), where the numbers in brackets denote the orienta- 
tion of the pure cation layers, A stands for Eu, and B for 
Gd. It is clear that no particular superlattice direction 
is preferred. 

For the region above x = 0.5, our phase diagram pre- 
dicts phase separation into the Lli phase and pure GdO. 
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TABLE III: Formation energies and the energetic depths of the unambiguously established (regular font) and tentatively 
predicted (italic) ground-state structures for Eui-^Gd^O and Eui-^Gd^S systems. A space group notation or a Strukturbencht 
designation of the cation ordering (if available) is given in parenthesis after the unit cell formula of each compound. AJ/gga, 
Ai/flxc, and AHce are the formation enthalpies obtained, respectively, from scalar-relativistic GGA+U with full relaxation, 
from GGA+U with a restricted relaxation in which the cations are fixed at the undistorted fee lattice with the Vegard-law 
lattice parameter (while the anions are allowed to relax), and from the CE fit. <5gga and 5ce are the energetic depths from the 
full relaxation and from the CE fit, respectively. AHqq A are similar to AHqga except the spin-orbit coupling was included. 
All energetic quantities are given in meV/cation. 



Formula unit 
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AHcca (AHrnn) 


A //fi xc 
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^GGA 


<^CE 






Eui x Gd x O: 










EuO (Bl) 





0.0 


0.0 


-10.6 






EuuGdOis (C2/m) 


1/15 


-33.0 


-32.8 


-31.2 


l.i 


0.5 


EuxiGdOvi (PI) 


1/13 


-36.8 


-36.5 


-33.9 


0.7 


0.0 


EugGdOg (Pilm) 


1/9 


-46.5 


-46.4 


-42.6 


3.2 


0.4 


Eu w Gd 2 0i2 (C2/c) 


1/6 


-54.0 


-53.8 


-55.7 


0.5 


2.7 


Eu 8 Gd 2 O w (PI) 


1/5 


-57.8 


-56.8 


-59.2 


1.0 


1.3 


Eu fi Gd 2 O s (C2/c) 


1/4 


-60.7 


-58.3 


-61.3 


0.4 


0.4 


Eu 8 Gd40i2 (C2/m) 


1/3 


-64.7 (-62.8) 


-61.3 


-63.7 


3.8 


2.9 


EuGd0 2 (Lli) 


1/2 


-61.5 (-59.9) 


-60.6 


-60.1 


12.9 


11.6 


GdO (Bl) 


1 


0.0 


0.0 


-2.6 










Eui_»Gd (r S: 










EuS (Bl) 





0.0 


0.0 


3.4 






EuGd 2 S 3 (C6) 


2/3 


-0.4 (-2.0) 


26.8 


-0.3 


0.4 


3.8 


GdS (Bl) 





0.0 


0.0 


3.6 







Since GdO, as has been mentioned above, has not been 
observed in the rocksalt structure considered here, we 
have further investigated its stability. In the zincblende 
structure its energy is found to be 56 meV/cation higher 
compared to rocksalt, but the equilibrium lattice param- 
eter 5.31 A is close to that (5.24 ±0.05 A) reported^ for 
zincblende GdO. The fact that experimental data sug- 
gest zincblende structure, while our calculation predicts 
it to be significantly less stable than rocksalt, may be 
due to the likely off-stoichiometry in experiment. Fur- 
ther, we considered the stability of GdO with respect to 
decomposition into Gd203 and metallic Gd. Gd20s oc- 
curs in three different structures^, cubic (80-atom unit 
cell) under 777 K, monoclinic (30-atom unit cell) from 
777 to 1727 K, and hexagonal above 1727 K. We used 
fully relaxed lattice constants of Gd203, Gd and GdO. 
Among the three phases of Gd2 03, the cubic phase is 
the most stable, its energy per formula unit being lower 
by approximately 90 meV/cation compared to the mon- 
oclinic phase. Rocksalt GdO energy was found to be 21 
meV/cation (per Gd atom) lower compared to the mix- 
ture of Gd203 and Gd. From this we can deduce that 
if the system is in equilibrium with an oxygen reservoir, 
GdO can only form in a narrow range of the oxygen chem- 
ical potential /i<> Indeed, using the following relations: 

Eodo ~ E Gd < fio < E Gd2 o 3 - 2E Gd o 
E Gd + E Gd2 o 3 = 3E Gd0 + 63 meV 

where the inequalities follow from the stability of GdO 
with respect to both reduction to Gd and oxidation to 
Gd2 03, we find that the double inequality is satisfied 
only in the 63 me V- wide range of /iq ■ In reality this range 



depends on temperature and is subject to the uncertainty 
in the calculated reaction enthalpy, but nevertheless this 
feature is consistent with the difficulty in stabilizing GdO 
experimentally. 

The calculated phase diagram for the oxide system is 
shown in Fig. [SJi. It is seen from Table IIIII that most 
of the ground states with x < 0.3 have rather small en- 
ergetic depths 5, suggesting that these orderings would 
only occur at very low temperatures. Indeed, our sim- 
ulations show that many of these phases only appear 
well below room temperature, so that the correspond- 
ing phase transformations are kinetically inaccessible. In 
fact, phase transformations occurring above T ~ 400 
K involve only the phases with x = 1/3 and x = 1/2. 
Therefore, for x < 0.3 we have only determined the ap- 
proximate boundary (i. e. the solubility limit) of the dis- 
ordered (Gd,Eu)0 phase, which is shown by a dot-dashed 
line in Fig. [5^. 

There are several interesting features in this phase di- 
agram. (1) A broad miscibility gap exists in the Gd-rich 
region with the critical point close to 1200 K and x « 0.7. 
(2) A continuous order-disorder transition occurs for the 
Lli phase, whose line terminates at a tricritical point 
(T w 600 K, x ~ 0.33) on the Eu-rich end and at a criti- 
cal endpoint at the Gd-rich end (T » 850 K, x « 0.51). 
Thermodynamics mandates that the slope of the solubil- 
ity line at the tricritical point should be different from 
the slope of the ordering line, but this difference is too 
small to be revealed in Monte Carlo simulations. On the 
other hand, the slopes of the binodals do not change at 
the critical endpoint, but their curvatures do. (3) There 
is a eutectoid triple point at T » 420 K at which the 
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disordered phase decomposes in a mixture of C2/m and 
Lli phases. (4) A point of equal concentrations caps 
the single-phase C2/m region; this point is close to the 
eutectoid point. 

Apart from the variations due to different phase transi- 
tions, the solubility of GdO in EuO grows approximately 
linear with temperature up to the critical point near 1200 
K, with a slope of about 0.06%/K. 



GdO 




GdS 



FIG. 5: Calculated phase diagrams obtained using Monte 
Carlo simulations from the cluster expansions for (a) 
Eu^Gdi-ajO and (b) EuxGdi-^S alloys. The phase labels 
indicate the ordering within the cation sublattice. In (a), the 
phases in the region bounded by the dot-dashed line have 
not been identified. In (b), the solid lines correspond to the 
paramagnetic high-temperature phases, while the dotted lines 
show the prediction of the T = ferromagnetic cluster expan- 



As discussed in Appendix [Cj the ordering transition 
for the prominent Lli ground state occurs at T or d = 840 
K and is second order around the stoichiometric EuGdC>2 
composition. The physics behind this ordering transition 
can be further illustrated by considering the CE Hamil- 
tonian [Eq. ([T]>] as a generalized Ising model, in which Eu 
atoms are represented by pseudo-spin 5=1/2 and Gd's 
by S = —1/2, and the pseudo-spins form an fee lattice. 
In the Lli structure, the nearest neighbor interaction is 
fully frustrated, while the second nearest neighbor inter- 
action is not frustrated (the second-nearest cation neigh- 



bors are always of dislike type). Thus, we can regard this 
structure as being formed by four interpenetrating sim- 
ple cubic lattices, all of which have AFM ordered pseudo- 
spins, and which are coupled only through longer-range 
interactions and through the order-from-disorder mech- 
anism. Indeed, the transition temperature in our MC 
simulations (about T or d = 840 K) is close to that of the 
simple cubic lattice under the assumption that its first 
nearest-neighbor interaction is equal to the second near- 
est neighbor interaction of the original fee lattice (920 K 
using the best known estimate of T c in the AFM Ising 
model from Ref. K59) . 



B. Eui-sGd^S 

The formation enthalpy diagram for the sulfide system 
(all compounds including GdS are assumed to be FM) is 
shown in Fig. [4]b. This sulfide system has a very different 
ground-state convex hull, compared to the oxide system. 
There is only one marginally stable compound at a; = 2/3 
with the formation enthalpy of only —0.4 meV/cation 
(see Table HHJ); this number is, in fact, smaller than the 
precision of our GGA+U calculations. This compound 
has the hexagonal C6 structure (sometimes referenced as 
"a2" in CE investigations); its structural parameters are 
given in Appendix |B] This structure is a superlattice 
composed of pure GdO and EuO (111) layers alternat- 
ing in a 2:1 pattern (similar in this respect to the Lli 
structure which has a 1:1 pattern). 

The computed phase diagram is depicted in Fig. [SJd. 
Dotted lines correspond to the FM CE, and solid lines 
to the corrected CE, designed to represent the PM phase 
equilibrium. This corrected CE was obtained by chang- 
ing the nearest-neighbor pair ECI in our FM CE from 
—6.46 meV to —7.37 meV, which results in the correct 
paramagnetic formation enthalpy of C6 EuGd2S3 (see 
Section IIVI and Appendix [XJ ■ 

In the FM phase diagram, there are two major fea- 
tures: a wide miscibility gap with a critical point at 
T w 810 K and x w 0.5, and a peritectoid triple point 
at T r; 560 K, at which the two disordered phases are 
in equilibrium with a new predicted C6 phase. The C6 
phase forms a very narrow single-phase region with con- 
centration slightly decreasing at elevated temperatures. 
The very small energetic depth S of the C6 phase is sig- 
nificantly overestimated by the CE fit (see Table IIII[) : 
we therefore expect that the temperature of the peritec- 
toid point is also overestimated. In fact, in the paramag- 
netic phase diagram (solid lines) the C6 phase is unsta- 
ble (5 < 0). Note, however, that the GGA+U value of 
5 (positive in the FM case and negative in the paramag- 
netic case) is comparable to the precision of our GGA+U 
calculations; we thus conclude that our accuracy is not 
sufficient to confidently select among the two scenarios 
shown in Fig [3b by the solid and dotted lines. 

We have also calculated the enthalpy of the reaction 
similar to the second line of Eq. [2] with the a-phase of 
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GCI2S3 and rocksalt GdS, which comes out at 2.8 eV. 
Thus, in equilibrium with a sulfur reservoir the rocksalt 
GdS is stable against both reduction and oxidation in 
the 2.8 eV-wide range of the S chemical potential. This 
is consistent with the fact that this phase is readily ob- 
tained experimentally. 

The overall shape of the sulfide phase diagram is much 
more symmetric compared to the oxide system, which is 
consistent with weaker triplet ECI's (see Fig. [3]). How- 
ever, the solubility limit of GdS in EuS is still greater 
than that of EuS in GdS, similar to the oxide system. 
The solubility limit of GdS in EuS is about 7% at 300 K 
and increases almost linearly at a rate of about 0.07%/K. 



VII. ANTON-MEDIATED DEFORMATIONAL 
INTERACTION 

The ground states with Lli or a closely related D4 
cation ordering (the latter having the same number of 
like and dislike neighbors as the Lli for any given neigh- 
bor distance), appear quite common for rocksalt chalco- 
genides: in addition to EuGd02, such ground states were 
recently predicted for a number of rocksalt telluridesj 60 i 61 
In the case of tellurides, D4 appears a more typical 
ground state, although Lli typically differs by only a 
few meV/cation in energy. Conversely, in EuGd02, D4 
is only 2.6 meV/cation higher in energy than Lli (which 
is sufficient for Lli to prevail up to T or( j = 840 K, due 
to the near identity of the entropic contributions to the 
free energy of both structures.) The appearance of Lli 
and D4 ground states has been related to the elastic soft- 
ness to a deformation along the [111] direction^ 2 which 
is indeed typical for rocksalt compounds, generally hard- 
est along the [100] cation-anion bond direction. Further, 
it appears that the other high-T or( j phase that we have 
discussed, namely the tentatively predicted C6 EuGd 2 S3 
phase, follows the same elastic trend, since it is a (111) 
superlattice just like Lli. 

It may not be immediately clear, however, why Lli 
EuGdS2 is not stable. More generally, while similar 
electronic structures for both oxides and sulfides are ex- 
pected from the isovalent electronic configurations (and 
confirmed by first-principles calculations), their ground- 
state formation enthalpy diagrams are very different: the 
oxides exhibit a convex hull with large and negative for- 
mation enthalpies, but the sulfides have positive forma- 
tion enthalpies for all ordered compounds except one at 
x = 2/3 which is close to zero (Fig. 2] and Table 1111]) . We 
shall now demonstrate that it is possible to understand 
the origin of this difference, as well as to get a deeper 
understanding of the elastic mechanism leading to the 
predominance of Lli and D4 in rocksalt chalcogenide al- 
loys, by examining the mechanism of atomic relaxation 
in these alloys. 

Let us decompose the formation enthalpy in two parts: 
the "unrelaxed chemical" part AH c h cm which is calcu- 
lated for the undistorted (cubic) lattice at the lattice pa- 



TABLE IV: Decomposition of the formation enthalpy for Llo 
and Lli structures in the unrelaxed chemical (Ai/ c h om ) and 
relaxation (AH re i) contributions (see text). Afl^J" an ^ 
AHord are the corresponding differences between the Llo and 
Lli structures. All enthalpies are given in meV/cation. 



Composition 


Structure 


^•^"chem 


AH re i 


A rjchcm 




EuGd0 2 


Llo 
Lli 


63.1 
37.1 


-2.4 
-96.5 


26.0 


94.1 


EuGdS 2 


Llo 
Lli 


88.1 
62.6 


-8.0 
-53.0 


25.5 


45.0 



rameter avL given by the Vegard law, and AH re i due to 
the additional full relaxation from these ideal Vegard po- 
sitions (thus, AH = AiJchem + AH re i). This decomposi- 
tion is shown in Table ITVI for the Ll and Lli structures, 
which straddle almost the entire range of formation en- 
thalpies at x = 0.5 for both oxide and sulfide systems. 
The difference AH or d between the enthalpies of the Llo 
and Lli structures gives the characteristic ordering en- 
thalpy; its decomposition in the unrelaxed chemical and 
relaxation parts is also included in Table IIV1 

First, we see that the "chemical" part of the ordering 
energy Aff™j™ is almost identical for both oxide and 
sulfide compounds, reflecting the fact that they are iso- 
electronic and their bonding properties are therefore very 
similar. Second, the values of Ai/chem arc significantly 
larger for both sulfide compounds compared to the oxide 
ones. This feature can be explained by a notably larger 
lattice mismatch in the sulfide system (see TablelJ), which 
forces the Gd-S and Eu-S bonds to deviate further away 
from their equilibrium lengths. 

For the tetragonal Llo cation ordering, symmetry for- 
bids the relaxation of all internal coordinates; only the 
lattice parameters are allowed to relax. It is seen from 
Table QV] that this relaxation yields only a relatively small 
energy gain. However, in the rhombohedral Lli struc- 
ture, the coordinates of the oxygen atoms are not fixed by 
symmetry. The relaxation of these coordinates makes the 
dominant contribution to AH re i, which is almost twice 
as large in the oxide system as it is in the sulfide sys- 
tem. Moreover, for Lli GdEu02 this internal relaxation 
overcomes the positive A-ffchcm contribution and makes 
the formation enthalpy large and negative. In the sulfide 
system AH re i is about two times smaller and does not 
fully overcome the Ai7 c h om term. 

Thus, we conclude that the Lli structure is strongly 
favored with respect to Llo because it allows the anion 
atoms to relax. Specifically, consider the octahedral cage 
occupied by an anion atom. In the Llo structure the Gd 
and Eu layers are stacked along the [001] direction; both 
inequivalent anion sites (in the Gd and Eu layers) are lo- 
cated at inversion centers. Thus, the opposite vertices of 
the octahedral cage are occupied by like cations (either 
both Eu or both Gd). On the other hand, in the Lli 
phase (or in the closely competing D4 phase) this inver- 
sion symmetry is fully broken in the sense that all the 
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opposing vertices of all octahedral cages are occupied by 
unlike cations. Since Gd and Eu ions have notably differ- 
ent radii, the anion atoms shift towards Gd. The struc- 
ture with the largest possible degree of inversion symme- 
try breaking maximizes the ability of the system to gain 
energy from anion relaxations. The opposite vertices of 
the octahedral cage are separated by the next-nearest 
neighbor distance in the cation sublattice. Thus, the re- 
laxation of the anion atom connecting the next-nearest 
neighbor pair is the mechanism generating the dominant 
positive ECI for this pair (see Fig. [3]). The fact that 
this ECI is smaller for sulfides is fully consistent with the 
data in Table ITVl and is not surprising because this anion- 
mediated interaction mechanism is sensitive to the anion 
size. Indeed, one can argue that the problem of finding 
the ground state with such interaction is largely equiv- 
alent to finding the most favorable packing of atoms to 
minimize the deviations from the optimal bond lengths. 

In order to understand the relative importance of the 
anion and cation relaxations, we have recalculated the 
formation enthalpies of all compounds by restricting the 
atomic relaxations to anions only and fixing the cations 
to the sites of an ideal fee lattice with the volume given 
by the Vegard law. These values are listed in Table IIIII 
as Ai?fi xc . For all oxide compounds the formation en- 
thalpies obtained in this way agree very well with the 
fully relaxed formation enthalpies (AH), indicating that 
anion relaxation in these structures are unrestricted by 
symmetry whereas the cation relaxations are insignifi- 
cant. For the EuGd2S3 compound in the sulfide system, 
however, the corresponding error is large. The reason is 
that this compound is a 2:1 layered superlattice, which 
makes it possible to adjust the interlayer Gd-S and Eu-S 
bond lengths by changing the cation layer separations. 
This is not possible in the Lli structure, which is a 1:1 
layered superlattice, because all the cation layer separa- 
tions are equal there. 



VIII. SUMMARY 

We have computed the temperature-phase diagram 
of two isovalent alloy systems Eui^^Gd^O, Eui-^Gd^S 
by using first-principles calculations combined with the 
standard cluster expansion and Monte-Carlo simulations. 
Very different ground-state convex hulls are obtained for 
the two systems: the oxides form ordered compounds 
with large and negative formation enthalpies, but sul- 
fides have only one marginally stable compound. The 
dominant configurational cation interaction comes from 
the second-nearest neighbor pair and is mediated by the 
relaxation of the O atom lying in between. The difference 
between oxides and sulfides is attributed to the difference 
in the anion size. Gd has a high equilibrium solubility 
in EuO and EuS even at room temperature, which in- 
dicates that rather heavy doping is possible without the 
precipitation of secondary phases. 
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Appendix A: Paramagnetic energies of sulfide 
compounds 

For pure EuS, GdS, and C6 EuGd2S3 structures pre- 
dicted to be stable by the FM CE, we investigated the 
paramagnetic phase regime by fitting magnetic configura- 
tional energies to the classical Heisenberg model. For the 
pure GdS, up to the third- nearest- neigbor exchange in- 
teractions were fitted using four magnetic configurations: 
(i) FM, (ii) AFM type I (AFMi) along the [001] direc- 
tion with alternating spins, (iii) AFM type II (AFMn) 
along [111] with alternating spins, and (iv) AFM type 
III (AFMin) along [001] with two layers of alternating 
spins. As expected, we found that AFMn is the ground 
state for GdS with the mean-field transition temperature 
of 63 K slightly above the experimental value 58 K. The 
fitting produces the paramagnetic formation enthalpy of 
—5.3 meV relative to the FM configuration. 

The paramagnetic formation enthalpy of the pure EuS 
was also obtained using the same magnetic configura- 
tions as for the pure GdS. The paramagnetic formation 
enthalpy turns out to be 1.3 meV above the FM phase 
with the mean-field value of the FM ordering tempera- 
ture 10.2 K, somewhat less than the experimental value 
of 19 K. 

For the C6 EuGd2S3 structure, only the nearest- 
neigbor interactions were considered but decomposed 
into four distinct types of Heisenberg exchange param- 
eters due to the layered structure and two cation species. 
In the [111] direction, the C6 structure establishes an 
A2B1 type superlattice, where each plane normal to the 
direction is composed of only one cation species. This 
superlattice structure allows two types of intra-layer ex- 
change interactions, Gd-Gd and Eu-Eu, and two types 
of inter-layer interactions, Gd-Gd and Gd-Eu. A to- 
tal of nine magnetic configurations were constructed in 
the 6-cation (1x1x2) supercell doubled along the [111] 
direction of the 3-cation unit cell, and also in the 6- 
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cation (2x1x1) supercell doubled in the normal direction 
to [111]. The magnetic configurations included, in the 
1x1x2 supercell, AfAfBtAtAtBt, AtA4.BtA4.AtB4., 
A4AtB4A4AtBt, A4A4BtAtAtB4, AtA4B4A4AtBt, 
AtAtBtA4AtB4, At AtB4 At AtB4 , and in the 2x1x1 
supercell, one of the two cations in either A or B layer 
has its spin flipped, while all other cation spins are kept 
opposite in directions. Its paramagnetic energy deter- 
mined by the fitting with the four exchange interaction 
parameters is slightly lower by 0.6 meV than the FM 
configuration. 

Appendix B: Predicted stable compounds 

The predicted ground-state structures with reasonably 
high ordering temperature are described below with their 
atomic positions fully relaxed using VASP. Additionally, 
for each structure, its crystallographic space group is 
identified for clarification. 



1. Eu g Gd 4 Oi2 

Space group #12 : C2/m, base-centered monoclinic 
Primitive unit-cell: a=12.507 A, b=7.187 A, c=8.766 A, 
a=90.00°, £=90.24°, 7 =90.00° 
Wyckoff positions: Eu(2i)=(0.084,l/2,0.339), 
Eu(2i)=(0. 7503,1/2,0.995), Eu(4j)=(0.832,0. 750,0.330), 
Gd(2g) = (0,0. 251,0), Gd(2i)=(0.417,l/2,0.667), 
O(2i)=(0.075,l/2,0.851), O(4j)=(0. 153,0. 237.0.168), 
O(2h)=(0,0. 224,1/2), O(2i)=(0. 260,1/2,0.512), 
O(2i)=(0.432, 1/2,0.159). 

2. EuGd0 2 

Space group #166: R3m, trigonal (rhombohedral) 
strukturbericht designation of cation order: Lli 
Primitive unit-cell: a=3.536 A, b=3.536 A, c=17.741 A, 
a=90.00°, £=90.00°, 7 =120.00° 
Wyckoff positions: Eu(lb)=(0,0,l/2), 
Gd(la) = (0,0,0), O(2c) = (0,0,0.260). 



Primitive unit-cell: a=4.041 A, b=4.041 A, c=9.959 A, 
a=90.00°, £=90.00°, 7 =120.00° 
Wyckoff positions: Eu(lb) = (0,0,l/2), 
Gd(2d) = (l/3,2/3,0.155), S(la)=(0,0,0), 
S(2d)=(l/3,2/3,0.681). 

Appendix C: Order-disorder transition in Lli 



In finite-size MC simulations there is no formal dis- 
tinction between a first- and second-order transition, and 
the nature of the transition can be unambiguously deter- 
mined only from a finite-size scaling analysis. Although 
we did not perform such an analysis, strong evidence in 
favor of the second-order character of the ordering tran- 
sition in Lli EuGd02 is revealed by the behavior of the 
enthalpy and the heat capacity. A MC heating simula- 
tion was performed starting with the Lli structure at 
T = K. Figure [5] shows the temperature dependence of 
the heat capacity, as well as the total enthalpy, along the 
composition line corresponding to a constant chemical 
potential /1 = 50 meV (the concentration is shown in the 
inset). The Lli phase persists up to 840 K, where a con- 
tinuous order-disorder transition is indicated by the char- 
acteristic peak of the heat capacity, as well as the contin- 
uous change of the enthalpy. Away from the Euo.5Gdo.5O 
composition, the Lli ordered and the disordered phases 
maintain the same features, although the heat capacity 
peak is reduced. 



0.003 




Temperature (K) 



3. EuGd 2 S 3 

Space group #164: P3ml, hexagonal 
strukturbericht designation of cation order: C6 



FIG. 6: (Color online) Monte Carlo heating simulation of the 
Lli-to-disordered phase transition for Euo.5Gdo.5O. Forma- 
tion enthalpy and heat capacity are shown by black and red 
lines, respectively. The inset shows the temperature depen- 
dence of the composition x. 
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